function [A, b] = GetAffineMatrixb(IDx, IDy, Ix, Iy)

%%
%initialize affine transform parameters
A = zeros(2,2);
b = zeros(2,1);
Ab = zeros(6,1);

%ID and Ixy are (n x 2) matrix
IDx1d = IDx(:);
ID = zeros(size(IDx1d,1),2);
ID(:,1) = IDx(:); ID(:,2) = IDy(:);

Ix1d = Ix(:);
Ixy = zeros(size(Ix1d,1),2);
Ixy(:,1) = Ix(:); Ixy(:,2) = Iy(:);

%SIDx = sum(IDx(:));
%SIDy = sum(IDy(:));

SIx = sum(Ix(:));
SIy = sum(Iy(:));

%testmax = zeros(size(IDx,1),size(IDx,2))+128;
%testmin = zeros(size(IDx,1),size(IDx,2))-127;
%OverlapMatrix = (IDx < testmax) .* (IDx > testmin);
%will get the right translation 
Sigma = min(size(IDx,1),size(Ix,1))*min(size(IDx,2),size(Ix,2));
%Sigma = sum(OverlapMatrix(:));

%%
Q11 = [1 0; 0 0]*ID';
Q12 = [0 0; 1 0]*ID';
Q21 = [0 1; 0 0]*ID';
Q22 = [0 0; 0 1]*ID';

%width = size(Ix,2);
%height = size(Ix,1);
%Op = sum((IDx>0).*(IDx<(width-1)).*(IDy>0).*(IDy<(height-1)));

%initialize parameter matrix
Am = zeros(2,2);
Bm = zeros(2,1);

I1Q11 = ID(:,1)'*Q11';
I2Q11 = ID(:,2)'*Q11';
I1Q12 = ID(:,1)'*Q12';
I2Q12 = ID(:,2)'*Q12';
I1Q21 = ID(:,1)'*Q21';
I2Q21 = ID(:,2)'*Q21';
I1Q22 = ID(:,1)'*Q22';
I2Q22 = ID(:,2)'*Q22';

% Am(1,:) = [Sigma 0 sum(ID(:,1)) sum(ID(:,2)) 0 0];
% Am(2,:) = [0 Sigma 0 0 sum(ID(:,1)) sum(ID(:,2))];
% Am(3,:) = [sum(Q11(1,:)) sum(Q11(2,:)) I1Q11(1) I2Q11(1) I1Q11(2) I2Q11(2)];
% Am(4,:) = [sum(Q12(1,:)) sum(Q12(2,:)) I1Q12(1) I2Q12(1) I1Q12(2) I2Q12(2)];
% Am(5,:) = [sum(Q21(1,:)) sum(Q21(2,:)) I1Q21(1) I2Q21(1) I1Q21(2) I2Q21(2)];
% Am(6,:) = [sum(Q22(1,:)) sum(Q22(2,:)) I1Q22(1) I2Q22(1) I1Q22(2) I2Q22(2)];
% 
% Bm(1) = SIx;
% Bm(2) = SIy;
% Bm(3) = sum(sum(Q11'.*Ixy));
% Bm(4) = sum(sum(Q12'.*Ixy));
% Bm(5) = sum(sum(Q21'.*Ixy));
% Bm(6) = sum(sum(Q22'.*Ixy));

%this set of parameters work
Am(1,:) = [sum(Q11(1,:)) sum(Q11(2,:))];
Am(2,:) = [sum(Q22(1,:)) sum(Q22(2,:))];
%the bug is '(sum(sum(Q22'.*Ixy))'
Bm(1) = (sum(sum(Q11'.*Ixy)) - I1Q11(1) - I2Q11(2)); 
Bm(2) = (sum(sum(Q22'.*Ixy)) - I1Q22(1) - I2Q22(2));

%%

%%
%the first two equations
%b(1) = (SIx - sum(ID(:,1)))/Sigma;
%b(2) = (SIy - sum(ID(:,2)))/Sigma;
%%
%the second two equations
%b(1) = (sum(sum(Q11*Ixy)) - sum(ID(:,1)))/Sigma;
%b(2) = (sum(sum(Q12*Ixy)) - sum(ID(:,2)))/Sigma;
%%
%Solve linear system by direct matrix inversion
%Aminv = inv(Am);
Ab = (pinv(Am)*Bm)';
%%
%Solve linear system by direct matrix inversion
%Ab = (Am)\Bm;
% Ab = Bm'/Am;
b = [Ab(1) Ab(2)];
%%
%Solve linear system by matlab
%b = linsolve(Am,Bm);
A = [1 0; 0 1];